rm(list = ls(all = TRUE)) #limpiar workspace
setwd("~/Dropbox/ENVIPE/SIMULA_SECUESTROS")
setwd("C:/Users/martin.romero/Dropbox/ENVIPE/SIMULA_SECUESTROS")

#-----------------------------------------------------------------------
#Generar una poblacion de <pob> habitantes, que viven en hogares de dos
#-----------------------------------------------------------------------
pob <- 10000000 #Tamano de la poblacion
prop <- c(.9995,.0001,.0001,.0001,.0001,.0001) #proporcion de hogares en cada categoria
data <- matrix(c(0,0,  #Hogar con ninguna victima
                 1,0,  #Una victima
                 0,1,  #Una victima
                 1,1,  #Dos victimas
                 2,1,  #Tres secuestros
                 1,2), #Tres secuestros
               ncol = 2, byrow=TRUE)
data <- data[rep(1:6, (pob/2)*prop),]

# Cuantos secuestros hubo?
real <- sum(data)
real

#---------------------------------------------------------------------
#Supongamos que deseamos hacer una encuesta de <size> hogares
#---------------------------------------------------------------------
size <- 500 #Numero de hogares en nuestra muestra
sims <- 50000 #Numero de simulaciones

##Muestra de hogares
set.seed(123456)
samples <- replicate(sims, data[sample.int(nrow(data), size=size),], simplify=FALSE)

#Estimadores
select <- matrix(rep(c(1,0, 0,1),size/2), ncol = 2, byrow=TRUE) #Seleccionar al respondente
ponde.indiv <- (nrow(data)/size)*2
ponde.hogar <- nrow(data)/size
indiv <- unlist(lapply(samples, function(x) sum(x*select))) *ponde.indiv
hogar <- unlist(lapply(samples, function(x) sum(x))) *ponde.hogar

mean(indiv)
mean(hogar)

sum(indiv>real)/sims
sum(hogar>real)/sims

quantile(indiv, c(0.025, 0.975))
quantile(hogar, c(0.025, 0.975))

#---------------------------------------------------------------------
#Introducir error de medición
#---------------------------------------------------------------------
gen.error <- function(data, falso.positivo, falso.negativo) {
  #Convertir a vector (desagregar individuos)
  x <- as.vector(data)
  #Falsos positivos
  x[x==0][sample(sum(x==0), size=round(sum(x==0)*falso.positivo))] <- 1
  #Falsos negativos
  x[x> 0][sample(sum(x> 0), size=round(sum(x> 0)*falso.negativo))] <- 0
  #Convertir a matriz (agrupar de nuevo en hogares)
  out <- matrix(x, ncol = 2, byrow=TRUE)
  return(out)
}
samples.w.error <- lapply(samples, gen.error, 
                          falso.positivo=0.02, 
                          falso.negativo=0.50)
#Estimadores
indiv.error <- unlist(lapply(samples.w.error, function(x) sum(x*select))) *ponde.indiv
hogar.error <- unlist(lapply(samples.w.error, function(x) sum(x))) *ponde.hogar

mean(indiv.error)
mean(hogar.error)

quantile(indiv.error, c(0.025, 0.975))
quantile(hogar.error, c(0.025, 0.975))
#---------------------
#Graphs
#---------------------
mi.plot <- function(dist, main, xlim, real) {
  hist(dist, 
       main=main,
       xlab="Miles de secuestros", 
       xlim=xlim)
  abline(v=mean(dist), lwd=2, col="red")
}

pdf("Graph_Hogar_vs_Victima.pdf", 9,4)
par(mfrow=c(1,2))
xlim <- c(0,max(indiv, hogar))
mi.plot(indiv, main="Estimador individual", real, xlim=xlim)
mi.plot(hogar, main="Estimador por hogar", real, xlim=xlim)
dev.off()

pdf("Graph_Error.pdf", 9,4)
par(mfrow=c(1,2))
xlim <- c(0,max(indiv, hogar))
mi.plot(indiv.error, main="Estimador individual", real, xlim=xlim)
mi.plot(hogar.error, main="Estimador por hogar", real, xlim=xlim)
dev.off()